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d We extend the concept of optical flow to a dynamic non-Euclidean setting. 

1 i Optical flow is traditionally computed from a sequence of flat images. It 

is the purpose of this paper to introduce variational motion estimation for 
images that are denned on an evolving surface. Volumetric microscopy im- 
ages depicting a live zebrafish embryo serve as both biological motivation 
and test data. 
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co 1 Introduction 

t-H 

Advances in laser-scanning microscopy and fluorescent protein technology have 
•tH increased resolution of microscopy imaging up to a single cell level [11]. They 

y\ allow for four-dimensional (volumetric time-lapse) imaging of living organisms 

and shed light on cellular processes during early embryonic development. Un- 
derstanding cellular development often requires estimation and analysis of cell 
motion. However, the amount of data captured is tremendous and therefore 
manual analysis is not an option. 

The specific biological motivation for this work is to understand the motion 
and division behaviour of fluorescence-labelled endoderm cells in a zebrafish 
embryo. The marked cells are known to be located on the surface of the embryo's 
yolk sac, where they form a monolayer [17] . Loosely speaking, they only sit next 
to each other but not on top of each other. Moreover, the yolk sac deforms over 
time; see Fig. [T] 
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Figure 1: Sequence of embryonic zebrafish images. The curved mesh represents 
a section of the yolk sac's surface. Depicted are frames no. 30, 45, 55, and 60 
of the entire sequence. All dimensions are in micrometer (/im). See Sec. 4.1 for 
more details on the microscopy data. 



We take these biological facts into account and restrict our attention to the 
analysis of cell motion on the yolk's surface. With this approach it is possible 
to reduce the amount of data by one space dimension. The resulting problem 
consists in the estimation of motion of brightness patterns that are restricted 
to an itself moving surface. We approach this problem by adapting the classi- 
cal concept of optical flow to the present setting, where the image domain is 
both non-Euclidean and dynamic. Note that due to the monolayer structure 
cell occlusions cannot occur. This makes the optical flow field a more reliable 
approximation to the true motion field. 

Our contributions in the field of optical flow are as follows. First, we for- 
mulate the optical flow problem on an evolving two-dimensional manifold and 
give two equivalent ways of linearising the brightness constancy assumption 
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(Sees. 2.1 and 2.2). One uses a parametrisation of the evolving surface, the 



other one is parameter-independent. Second, we use a generalisation of the 



Horn-Schunck model to regularise the optical flow field (Sec. 2.3). For a given 
global parametrisation of the evolving surface, we solve the associated Euler- 
Lagrange equations in the parameter domain with a finite difference scheme 
(Sec. [3]). Finally, we apply this technique to obtain qualitative results from the 
afore- mentioned zebrafish data (Sec. [4]). Our experiments show that the optical 
flow is an appropriate tool for analysing these data. It is capable of estimating 
global trends as well as individual cell movements and, in particular, it is able 
to indicate cell division. 

Related work. Optical flow is the apparent motion in a sequence of images. Its 
estimation is a key problem in Computer Vision. Horn and Schunck [5] were the 
first to propose a variational approach assuming constant brightness of moving 
points and spatial smoothness of the velocity field. Since then a vast number of 
modifications has been developed. See [1] for a recent survey. 

Miura [13] observed that until 2005 optical flow has been mostly disregarded 
as a method for motion extraction in cell biological data. Since then, a few 
articles have explored this direction: Melani et al. [12] and Hubeny et al. [6] 
extended variational optical flow methods to volumetric images to obtain 3D 
displacement fields. In the former article, the resulting algorithm is also applied 
to zebrafish microscopy data. Quelhas et al. [14] use optical flow to detect cell 
divisions in a live plant root. However, they work with 2D (plus time) data only. 
Therefore, their approach suffers from errors caused by 3D off-plane motion. 

Clearly, certain natural scenarios are more accurately described by a velocity 
field on a non-flat surface rather than on a flat domain. With applications to 
robot vision, Imiya et al. |7J [16] considered optical flow for spherical images. 
In a more general setting, Lefevre and Baillet [10] extended the Horn-Schunck 
method to 2-Riemannian manifolds and showed well-posedness. They solve the 
numerical problem with finite elements on a surface triangulation. In all of the 
above works the underlying imaging surface is fixed over time, while in this 
paper it is not. 



2 Optical Flow on Evolving Surfaces 
2.1 Brightness Constancy 

Let M. t C M 3 , t G I = [0,T), be a compact smooth two-dimensional manifold 
evolving smoothly over time. We assume the velocity to be unknown. Moreover, 
denote by / a scalar time-dependent quantity defined on the surface 

/: U ( M * X {*» M ' 

tei 

We begin with a Lagrangian specification of the optical flow field. That is, 
for every starting point Xq G A4q we see k a trajectory where the data / are 
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conserved. More precisely, we want to find a function 

7: Mo x I -> [jMt, 

tei 

such that 

1. j(x ,t) e M t for all tei, 

2. j(-,t) is a diffeomorphism between Mo and Mt for all tei, 

3. 7(-,0) = ld Mo , 

is fulfilled and which satisfies a "brightness" constancy assumption (BCA) 

/(a? , 0) = f(nf(x ,t),t), for all (a? , G M x J. (1) 

In classical optical flow computations it is common practice to linearise the 
BCA by taking its time derivative and to solve the resulting equation for the 
Eulerian unknown 7^] We also take this route, but differentiation of / is more 
involved. Observe, for example, that for an arbitrary to e I and x e Mt the 
usual partial derivative 

d t f{x, t ) = lim - (f(x, t + h)- f(x, to)) 

is not well-defined, simply because, in general, x is not an element of Mt +h 
for all h ^ 0. 

In the next section we linearise in two different ways. First, we use a 
global parametrisation to pull the data back to a fixed reference domain and 
linearise afterwards. In our second approach we borrow some notions from 
continuum mechanics [2 to directly linearise ([!]). 

2.2 Linearisation 

Linearisation after pull-back. Let Q e M 2 be a compact domain and 

x: Q x / R 3 , (xi,x 2 ,t) = (x,t) \-> x(x,t) e M t 

be a parametrisation of the evolving surface. Denote by / the coordinate rep- 
resentation of /, that is, 

f(x,t)=~f(x(x,t),t) (2) 

and let 

13: n x / n 

be the coordinate counterpart of 7. This means, if we let xo = x(xo,0), then 
f3(xo,t) gives the coordinates of j(xo,t) e Mt in ft (see Fig.|2|. In other words, 
we have the identity 

7OO0, 0), t) = x(/3(xo, t),t), for all (x , t) eft x I. (3) 

1 To simplify expressions we use Newton's notation for those time derivatives that corre- 
spond to actual velocities, for example 7 = dtj. 
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Now, from 0, ^ and ([3) we get 

f(x ,0) = f(x ,0) 

= fa(xo,t),t) 

= f(x(P(x ,t),t),t) 

= f(P(x ,t),t), 

which is a coordinate version of the BCA. After differentiation with respect to 
t it becomes 

V 2 / • $ + d t f = 0, (4) 

where V 2 = (<9i,d2) T is the two-dimensional spatial gradient. Note that the 
last equation is nothing but the classical OFC for Euclidean data / and a dis- 
placement field $. 

Direct linearisation. We turn to our second derivation. While, as pointed 
out above, the partial derivative dtf is undefined in general, it does make sense 
to differentiate / following the surface movement. Let y be a point on Mt and 
£: t \-> £(£) G Mt an arbitrary smooth trajectory through the evolving surface 
satisfying £(£o) = V- Now we can compute 

^/(».*o) = lim ^ (/(C(*o + *0,*o + h)- f(y,t )) 

to obtain a valid derivative of /. Since this time derivative only depends on the 
vector v = £(£o)? we denote it by d^f. A natural candidate for a trajectory along 
which to differentiate is given by the parametrisation £(£) = x(x,t). Another 
possible choice would be a trajectory that is normal to Mt - The resulting 
normal time derivative is accordingly denoted by d™f. 

Finally, we also need the surface gradient V mJ If is a smooth extension 
of / to an open neighbourhood of y G Ait in then the surface gradient of F 
at y is defined as the projection of the three-dimensional spatial gradient V 3 F 
onto the tangent plane to Mt 

V^F = V 3 F - (V 3 F • n)n, 

where n is the unit normal to Mt - The surface gradient only depends on 
the values of F on the surface; see e.g. [4j p. 389]. Thus, V Mf — V 'mF 1S 
well-defined. 

The spatial and temporal derivatives of / introduced above are related in a 
simple way. As shown in [2], they satisfy the equality 

(5) 

= V^/'^tan + d^/, 

where x tan is the tangential surface velocity, that is, the projection of x onto the 
tangent plane to Mt - This decomposition of dff into normal and tangential 
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components is clearly valid for any trajectory in place of and therefore in 
particular for the unknown 7. This means we can use ^ in order to differentiate 
the BCA with respect to t. The resulting OFC reads 

V^/~-7tan + dr/~ = 0. (6) 

Discussion. We conclude this section with a brief comparison of the two OFCs 
derived above. We start by showing how to obtain Q from ^ and vice versa. 
To this end we again assume the existence of a global parametrisation and 
rewrite all quantities in ^ in terms of x. First observe that, by ([3]), the velocity 
of 7 equals the surface velocity x plus a purely tangential component 

7 = x + J/3, 

where J = (d\X &2x) is the Jacobian matrix of x with respect to x. On the 
other hand, by ([5]), the normal time derivative is equal to the time derivative of 
/ following x minus its tangential component 

d?f = dff-V M f-x. 

Using the last two equations to rewrite the left-hand side of ^ yields 

V M f ■ 7 + d?/ = V M f ■ (x + J/j) + dff - W M f ■ x 
= V M f-JP + dff, 

which is already the left-hand side of Q in terms of /. It only remains to observe 
that dff = d t f and to replace the surface gradient V ' m$ by its coordinate 
expression J^ _1 V 2 /, where g = J T J is the coefficient matrix of the Riemannian 
metric; see e.g. [9]. 

The OFC is a valid approximation to the BCA only in case of data with a 
sufficiently fine temporal resolution or, equivalently, for data with small displace- 
ments. This is a well-known limitation in classical optical flow computations. 
We highlight that for the optical flow problem on evolving surfaces there are 
two different limitations, depending on whether one considers Q or (J6|. In (4) 
the unknown is /3, while in ([6| it is 7 tan = #tan + J$- This means that (4) 
constrains the motion relative to the tangential surface velocity cb ta n, while (6) 
constrains the absolute tangential motion. 

Fortunately, the temporal resolution of our microscopy data is sufficiently 
high. We therefore use a natural global parametrisation (see Sec. [3]) to pull 
the data back to the Euclidean plane and solve Q. However, equation ^ is 
independent of any parametrisation. It can thus serve as a starting point for 
alternative numerical approaches. 

2 . 3 Regular isat ion 

From now on we fix an arbitrary to G / and turn to the actual solution of the 
parametrised OFC for (ui(x), U2(x)) T = u{x) = $(x,to). Recall that with this 
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o > n 



x(.,t) 



M > M t 



Figure 2: Commutative diagram describing the relation between unknowns f3 
and 7. 



notation u contains the coefficients of the tangential vector field u = J/3 with 
respect to the tangential basis (d\X,d2x) of Mt - Note also that, by fixing to, 
there is no more time-dependence in our problem which makes it effectively an 
optical flow problem on a static surface. Hence we omit any reference to to from 
now on and write M. instead of Ait - 

The sought vector field is underdetermined by the OFC alone. We over- 
come this by minimizing a functional that penalises violation of the OFC while 
imposing an additional smoothness restriction on u. More precisely, we use 
a recent extension of the original quadratic Horn-Schunck regularisation to a 
Riemannian setting [10]. Basically, they propose to minimise 

Ql 1 1 ~ •~ii2 In ii2 In ii2 

-u + d*f\\ L2{M) + 2\\Viu\\ L2(M) + ^2u\\ L2{My (7) 

Here, a > is the regularisation parameter and V ' jU denotes the covariant 
derivative of u along tangential basis vector djX, which is the conventional 
derivative projected onto the tangent plane of M. Using parametrisation x and 
the corresponding Christoffel symbols T % - k (see Sec. [3]) it can be written as 

2 22 
VjU = ^2 D ij d i x = ^2 {pi Ui + ^ TZ jk u k)diX, j = l,2. (8) 

i=l i=l k=l 

Additionally, for a potentially vector- valued function w on M, in Q we used 
the notation 



\M\lhm) 



= / \w\ 2 dM = / \w o x\ 2 y/detgdx. 
Jm Jn 



Rewriting ^ in terms of coordinate expressions, we arrive at the equivalent 
functional 

£(u) = l [ \a(V 2 f-u + d t f) 2 + \\D\\l]^^gdx, (9) 

where ||-D||f is the Frobenius norm of the covariant derivative matrix D = (Dij) 
as defined in (l8l). 
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3 Numerical Solution 



We solve the problem of minimising functional £ via its associated Euler-Lagrange 
equations. Regarding the integrand of £ as a function G(x, u, V 2, ui, V 2 ^), they 
read 

G Ul = d\Gd lUl + d^Gd^ux 
G U2 = diGo lU2 + d2Gd 2 u 2 i 

where subscripts of G denote partial derivatives. The resulting pair of linear 
PDEs is of the form 

Ai^i = V 2 u\ ■ c + V 2 U2 ■ d + u ■ bi + a\ 

AU2 = V 2 U 2 ' C + V 2 Ui • d J rU-b2 J r a 2 . 

The coefficient vectors a, 61,62, c, d are rather lengthy functions of the data / 
and metric tensor g, which is why we do not write them out in full here. Letting 
Vt = (0, l) 2 for simplicity, the natural boundary conditions of the variational 
problem are 

d j u i + y 52r) k u k = 0, for xj G {0, 1}, (11) 

k 

where z, j G {1,2}. In case of a flat manifold, e.g. M, = the Euler-Lagrange 
equations (10) reduce to those of the original Horn-Schunck functional and the 
boundary conditions become the usual homogeneous Neumann ones. For more 
details on the calculus of variations we refer to [3] . 

Due to the nature of the microscopy data (see Sec. 4.1 and Fig. [T]), the 
manifold A4 t modelling the deforming yolk sac is a surface with boundary that 
is most easily parametrised as the graph of a function z : ft x I — >> R. Hence, 
we set x(xi,X2,i) = (x\, X2, z(xi, X2, t)) T . Accordingly, for the metric we get 

g = h + V 2 zV 2 z T , det g = 1 + | V 2 z| 2 , 

where I2 G R 2x2 is the identity matrix. The Christoffel symbols turn out to be 

Ti jk = 1^2^ (djgkm + d k g mj - d m g jk ) = - 

Partial derivatives of z and of the projected data / were approximated by 
central differences. The system (10) with boundary conditions (11) was then 
solved with a standard finite difference scheme. In the following section numer- 
ical results are presented. 



4 Experiments 
4.1 Data 

As mentioned before, the biological motivation for this work are cellular im- 
age data of a zebrafish embryo. Endoderm cells labelled with green fluorescence 
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were recorded via confocal laser-scanning microscopy resulting in time-lapse vol- 
umetric (4D) images; see [11] for the imaging techniques. This type of image 
shows a high contrast at cell boundaries and a low signal-to-noise ratio in gen- 
eral. Our videos were obtained during the gastrula period, which is an early 
stage in the animal's developmental process and takes place approximately five 
to ten hours post fertilisation. In short, the fish forms on the surface of a 
spherical-shaped yolk sac; see e.g. [8] for many illustrations and detailed expla- 
nations. For the biological methods such as the fluorescence marker and the 
embryos used in this work we refer the interested reader to [15] . The important 
aspect about endodermal cells is that they are known to form a monolayer dur- 
ing gastrulation [17], meaning that the radial extent is only a single cell. This 
crucial fact allows for the straightforward extraction of a surface together with 
a two-dimensional image of the stained cells. Since only a cuboid region of ap- 
proximately 860 x 860 x 340 fim of the pole region is captured by the microscope, 
this surface can easily be parametrised; cf. Sec. [3] The spatial resolution of the 
Gaussian filtered images is 512 x 512 pixels and all intensities are given in the 
interval [0,1]. Our sequence contains 77 frames recorded in intervals of 240 s 
with clearly visible cellular movements and cell divisions. 

4.2 Numerical results 

In the following, we present qualitative results and demonstrate the feasibil- 
ity of our approach. For every subsequent pair of frames we minimised the 
functional ([9| as outlined in Sec. [3] We chose grid size as well as temporal 
displacement as ft = 1 and the regular isat ion parameter was set to a = 10. 
For demonstration purpose we make use of the standard flow colour-coding [1] , 
which maps (normalised) flow vectors to a colour space defined inside the unit 
circle. It is easy to see that the same colours are valid all over the manifold due 
to the parametrisation. 

As representative candidates for this discussion we chose the displacement 
field between frames 57 and 58 for the following reasons. First, the surface is 
distinctly developed. Second, a considerable number of cells is present in the 
image and third, the interval contains cell divisions. Fig. [3j left, shows the 
colour-coded tangential vector field and the colour space whereas Fig. [3] right, 
displays the same motion field as computed in the parameter space. A visual in- 
spection of the dataset shows that cells tend to move towards the embryo's body 
axis, which roughly runs along the main diagonal in Fig. [3] right. Clearly, the 
velocity field is sufficiently smooth and suggests this behaviour in an adequate 
manner on a large scale. The expected change in orientation along the body 
axis is well represented by the colour shift from orange-yellow below the main 
diagonal to purplish blue in the region above. On the contrary, the choice of the 
regularisation parameter ensures that individual movements are well preserved 
as can be observed from the image. 

Fig. [4] gives a detailed view of the section outlined by a (red) rectangle in 
Fig. [3] right. This section was chosen because it depicts a cell division. Fig. [4j 
left, and Fig. [4] right, display the frames before and after the event, respectively. 
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Figure 3: Optical flow field between frames 57 and 58 of the sequence. Colours 
indicate direction whereas darkness of a colour indicates the length of the vector. 
Note that the colour circle has been enlarged for better visibility. 

Moreover, in Fig. [4j left, the velocity field is shown. From the raw data we 
observed that when a cell actually splits, the two daughter cells drift apart in 
a 180 ° angle with respect to the mother cell. The displacement field clearly 
shows the anticipated pattern caused by the diverging daughter cells. In Fig. [3j 
right, the event is point up by two areas which are coloured mutually opposite 
with respect to the colour space. Our results suggest that cell division can be 
indicated reasonably well by our model. 

5 Conclusion 

Aiming at efficient motion analysis of 4D cellular microscopy data, we gener- 
alised the Horn-Schunck method to videos defined on evolving surfaces. The 
biological fact that the observed cells move along an itself deforming surface 
allows for motion estimation in 2D (plus time). In the course of this work, we 
presented two ways to linearise the brightness constancy assumption and showed 
that one could be obtained from the other and vice versa. The resulting opti- 
cal flow constraint was solved by means of quadratic regularisation and verified 
on the basis of the afore-mentioned data. Our qualitative results suggest that 
both global trends as well as individual movements including cell division are 
well shown in the surface velocity field. However, so far we only laid the basic 
groundwork in terms of a mathematical model. 
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Figure 4: Detailed view of a cell division occurring between frames 57 (left) 
and 58 (right). All vectors are scaled and only every fourth vector is shown. 
Intensities are interpolated for smooth illustration 
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